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ABSTRACT 

In the solar convection zone acoustic waves are scattered by turbulent sound 
speed fluctuations. In this paper the scattering of waves by convective cells is 
treated using Rytov's technique. Particular care is taken to include diffraction 
effects which are important especially for high-degree modes that are confined 
to the surface layers of the Sun. The scattering leads to damping of the waves 
and causes a phase shift. Damping manifests itself in the width of the spectral 
peak of p-mode eigenfrequencies. The contribution of scattering to the line 
widths is estimated and the sensitivity of the results on the assumed spectrum 
of the turbulence is studied. Finally the theoretical predictions are compared 
with recently measured line widths of high-degree modes. 
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1 INTRODUCTION 

In the solar convection zone thermal instabilities give rise to turbulent convection as the pre- 
dominant heat transport mechanism. Turbulence is one of the unsolved puzzles of modern 
physics, and the lack of a fundamental theory of turbulence requires a variety of phenomeno- 
logical models. This constitutes a serious limitation to the predictiveness of many astrophys- 
ical models and very often leads to a built-in uncertainty (Canuto & Christensen-Dalsgaard 
1996). Theoretical modelling of convection in the Sun is extremely difficult. Solar convection 
is highly turbulent and occurs on a vast range of scales. The material is compressible and 
magnetic fields can play an important role. Consequently, various parametrizations have 
been invented to describe the properties of the convection. At the same time vast efforts 
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are being undertaken to numerically simulate solar convection. These simulations are still 
at an early stage and their results are not yet conclusive. One of the prime motivations for 
studying solar and stellar structure in great detail is to improve and test these models of 
convection. 

Statistically significant differences between the observed and theoretically predicted eigen- 
frequencies of the Sun persist (Gough et al. 1996). The differences between computed adia- 
batic eigenfrequencies of the standard solar model and measured p-mode frequencies for any 
given order n increase with frequency when scaled with the mode inertia (with the model 
frequencies being bigger). Higher frequencies of modes of a given radial order correspond to 
greater values of the spherical degree /, and hence to modes that remain increasingly confined 
to the surface layers of the Sun. This suggests that the cause for the discrepancy between 
theory and observation lies in an inadequate modelling of the surface layer of the Sun. It is 
this surface layer where convection is important and constitutes the greatest uncertainty in 
the model. 

Helioseismic inferences are not reliable unless details of stellar structure and wave propa- 
gation within the acoustic cavity are correctly treated. It has been noted elsewhere (Canuto 
& Christensen-Dalsgaard 1996) that the errors in the computed frequencies caused by un- 
certainties in the equation of state (EOS) are large compared to the observational errors. 
The variation in the average wave propagation properties inside the Sun caused by turbulent 
convection and temperature fluctuations constitute yet another process that will affect the 
eigenfrequencies. This influence is difficult to assess and has not received as much attention 
as the errors introduced by uncertainties in the EOS or the opacity. Since helioseismic in- 
ferences are dependent on the detailed treatment of turbulence, the Sun could serve as a 
laboratory to study turbulence under conditions unattainable on Earth. 

Turbulence enters the theory of solar oscillations in two places: it changes the model 
of solar structure mainly through the turbulent pressure that provides additional support 
against gravity, and secondly it interacts with the waves themselves. The waves interact 
with the turbulence via two mechanism: The turbulent velocity advects the acoustic waves 
(Brown 1984; Swisdak & Zweibel 1998) and the temperature fluctuations cause fluctuations 
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Probing Solar Convection 3 
in the refractive index and thus in the local phase speed of the waves. 

The good quality of solar oscillation data does not only allow an accurate determination 
of the eigenfrequencies but also contains detailed information such as the width of the 
spectral peak. 

The observed line width of acoustic modes is caused by absorption and other non- 
adiabatic effects (e.g., Gough 1980; Goldreich & Kumar 1991; Balmforth 1992), as well 
as by scattering by turbulent velocity fluctuations in the solar convection zone. 

The scattering contribution to the line width has been studied by Goldreich & Murray 
(1994) in the framework of normal modes. Using a simple model of scattering of standing 
waves in a box and within the strict limits of geometrical acoustics, they derived a scatter- 
ing width for radial modes of ~ uM"^ /T{{n + 1), where uj is the angular frequency of the 
mode, M = u/cis the Mach number of the turbulent velocity and n the order of the mode. 
However, geometrical acoustics is a poor approximation for non-radial high-degree modes, 
because at the top of the acoustic cavity the correlation length of the perturbations becomes 
comparable to the wavelength of the waves. 

In this paper the scattering contribution to the line width is calculated, explicitly taking 
diffr active effects into account. Recently measured line widths of high-degree modes are 
compared with results of this calculation. Thus, one of the aims of this study is to investigate 
to what extent the structure of solar convection can be probed by seismic methods. For this 
purpose a formalism to describe the propagation of waves through random fluctuating media 
is introduced. For the sake of clarity I start by deriving the formalism for a homogeneous 
medium before generalising it to inhomogeneous media in Sec. 2.2. 

2 WAVE PROPAGATION THROUGH FLUCTUATING MEDIA 

Let us write down a wave equation with a scattering term on the rhs: 

L%1) = (V^ + k^)^ = ey(r)^, (1) 

where V is sometimes called the 'scattering potential'. For a homogeneous medium eV{r) = 
2/c^5c(r)/c(r), where 6c is the rms of the sound speed perturbation. As the scattering is 
assumed to be weak, e is a small number. Thus one can solve equation (|l|) perturbatively. 
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An approximate solution can be written in terms of a Born expansion (e.g. Nayfeh, 1973) 
by writing 



i;=Y: (2) 



m=0 



Substituting into equation (|I|) yields 



(V2 + p)(V;o + eiJ^ + ...) = ey(r)(V'o + e^-i + •••)• (3) 
Equating equal powers of e yields for m = 0: 

L^o = (4) 

and recursively for m > 

Liprn = V{r)^IJrn-l- (5) 

Equations and (|^) can be solved successively using Green's function 

= J rf3riK(ri)^„_i(ri)G(r;ri), (6) 
where the Green function, G, is defined by 



(V2 + fc2)G(r;ri) = 5^(r-ri). (7) 

For a homogeneous medium ^(r; ri) = exp(i/c|r — ri|)/47r|r — ri|. From equation (P) one can 
write 

7/>i(r) = J dh,V{r^)G{r- ri)^o(ri), (8) 

and successively 

^^{r) = I dh,... I dh^V{r,)...V{r^)G{v2; r,)...G{r; r„)V^o(ri). (9) 

An alternative and possibly more intuitive technique to calculate phase and amplitude 
fluctuations, known as Rytov's method, starts by recasting the wavefield into an exponential 
(see, e.g., Ishimaru 1978): 
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^(r) = e*W. (10) 

Hence one can write 

V2^ = V^[V^- + V^^], (11) 

and equation (|1|) becomes 

V^* + ■ + P - e\/ = 0, (12) 

which is a nonhnear first-oder differential equation known as Riccati's equation. Denoting 
\l/o := ^1/ in the absence of fluctuations, i.e. when = 0, and writing \f' = \f'o + \E'i, one 
obtains 

V^^i + 2V^o • V^i = -[V^i ■ V^i - eV]. (13) 

With the identity 

V2(^o^i) = (V^^o)^! + 2^oV^i + V'oV^^i, (14) 
one obtains the following inhomogeneous equation for ipo'^i 



(V2 + P)(^o^i) = [V^i ■ V^i - eV]^o, (15) 
which can be solved using Green's function, yielding 

^i(r) ~ -1- / G(r - r')[V^i ■ V^i - eV]M^')dh'. (16) 
To first order, \E'i = in the integrand, so that 



^i(r) ~ — ^ / eVG{r - r')^o(r')rf'r'. (17) 
Thus the first Rytov solution is given by 



V^(r) = V^o(r)e'^^^^^ (18) 
It has been shown that the first term in the Rytov expansion is superior to the first Born 
approximation (Keller 1969), and I shall employ Rytov's technique in the following analy- 
sis. The two methods discussed above, the Born series and Rytov's technique, are the most 
commonly used techniques dealing with scalar wave propagation through random media. 
More recently, Samelsohn & Mazar (1996) have treated this problem using a path-integral 
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analysis on the basis of the parabohc wave equation. For a detailed review of stochastic wave 
propagation and scattering in random media, see, e.g. Ishimaru (1978) and Klyatskin (1980). 



Here it will be assumed that the sound speed fluctuations, and consequently the 
fluctuations are Gaussian with zero mean, i.e. (\E'i) = 0, where the brackets denote the 
mean over time. The mean of equation (plSj) has the form of a characteristic function and 
one can write 



where {"^D is the correlation function defined in equation (24) (see Munk & Zachariasen 
1976, Panchev 1971). 

Whenever waves interact with obstacles or inhomogeneities in the medium, diffraction 
effects may occur. Most problems in diffraction theory cannot be solved exactly, so that a 
number of techniques have been invented to find approximate solutions. One of these tech- 
niques is Huygen's principle, which states that every point of a wave front can be considered 
as a source of secondary wavelets which mutually interfere. The application of this princi- 
ple to diffraction problems leads to the Fresnel zone construction. The textbook example 
for the application of Fresnel zones is the problem of diffraction by a small aperture in 
an infinite screen as sketched in Fig. 1. A wave propagates from a point source at point S 
through an opening in an opaque screen to point R. According to Huygen's principle the 
total disturbance at R due to a source at point S is given by the integral over the aperture 



where s = SP and r = PR. In equation (^) I have neglected a so-called obliquity factor, 
which describes the angular variation of the secondary wavelets. Rays that connect S and 
R via a point P in the plane now have a phase at R that is different from the phase of the 
stationary ray. 

The envelopes of all the rays of length a + b + X/2, a + b + X, a + b + jX/2, divide 
the aperture into a number of zones, which are called Fresnel zones (a and b are defined in 
Fig. 1, A is the wavelength and j an integer). The contributions from successive Fresnel zones 
to the total wavefield alternate in sign. Those rays whose phase differs from the stationary 
phase by tt form the boundary of the first Fresnel zone in the plane of the obstacle. In our 




(19) 




(20) 
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textbook example the first Fresnel zone is a concentric circle centred on the optical axis. For 
a more detailed treatment see, e.g., Born & Wolf (1970). 

Generally, the first Fresnel zone in a plane intersecting the raypath is defined as the area 
bounded by all those rays whose phase differs by vr from the phase of the stationary ray which 
joins two fixed points. This can now be repeated for all planes intersecting the ray path and 
the entirety of these zones yields a tube surrounding the stationary ray path: this tube I call 
the Fresnel tube. The Fresnel zone provides an estimate of the size of an irregularity that 
would give rise to diffraction effects. If the size of the irregularity is much bigger than the 
Fresnel zone, diffraction effects are negligible. It is straightforward to calculate the Fresnel 
zone for waves in a homogeneous medium. Let me consider a wave propagating from S 
to R as before. If the wave now follows the path SPR instead of the stationary ray, the 
corresponding change in phase, 0, is given by 

50 = fc[(a2 + 5^2)1/2 ^ ^ ^^2]i/2 _ ^ ^ '<^±^sz^ (21) 
where k denotes 2tt/X. Hence the size of the Fresnel zone is given by 

2.1 Homogeneous medium 

I begin by considering the case of wave propagation in a fluctuating homogeneous medium. 
In this case the ray paths are straight lines and the turbulent correlation function is a 
function of separation only, i.e. 

e2(V(ri)V(r2))=C(ri-r2). (23) 

As stated in equations (|1^) and (jl^) the time-averaged wavefield at r due to a source at the 
origin is given by the unperturbed wavefield times exp(i(\E'f )), where 

2 

i^l) = (IQnY' I dh, [ dh2- ^— rC(ri - r^) 

J J ri|r — ri|r2|r — 

X exp[iA;(ri + |r — ri| — r)] exp[iA;(r2 + |r — r2| — r)]. (24) 
It is convenient to convert to relative and centre-of-mass coordinates: 



27[ab 
k{a + b) 



(22) 
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Ti - ra 



(25) 



and 



1 



R= - ri 



rs). (26) 

I assume that ( is negligible when f is bigger than the correlation length of the sound speed 
perturbations, which are taken to be smaller than R. Then one can expand the integrand in 
(24) in r/R yielding 



1 



167r2 



exp[2ik{R + |r 



i?2|r-R|2 

^^-^ 



ik 
T 



-2 



R| -r)] 
- (? ■ R)2 



R)]^ 



R ^-^rrll '''' 

where the hat denotes a unit vector. Assuming that the radius of the Fresnel zone is smaller 
than s, i.e. (r — s)/krs <^ 1, one can expand equation (27) in Rj_/R\\ about its stationary 
path, which is the straight line from to r. Here s stands for R\\. The subscripts _L and || 
denote the vector components perpendicular and parallel to r. Thus one can write 



Moreover, 



and 



i?+ |r - RI - r 



rV 



f2 - ff . R)2 p2 



R 



2s(r — s) 



(28) 



O 



(29) 



Substituting 
finds that 



- [r • (f - R)]- 
|r-R| 



O 



(30) 



into (27) and evaluating the integral over R_l by stationary phase, one 



d^n— — rexp[2iA;(/?+ |r-R| -r)l ~ — / 

Rp k Jo s[r — 



ds. 



(31) 



i?2|r — - - ■ ■ - ft Jo s[r — s) 

Neglecting terms of order r^R^/s"^ in equation (27) one is left with the integral over r which 

after substitution of (^) and (|30|) becomes 
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- TtV /'^^ / C(g) . s X exp 
Iott/c jo j s(r — s) 



ik (\ 

^ - + - 

4 Vs r 



(32) 



It is convenient to use the Fourier transform of the correlation function, which is defined as 



C(q) = / C(f)e'''-' d:^v. 
Substituting the inverse tranform into equation (|32D , using the relation 



(33) 



1 C(r||,0)df|| = (Svr)-^! C(0,q^)d2q^, 



and integrating over (i^r^ by stationary phase, one obtains: 



(34) 



(i^q^ C(0, q±) / ds exp 

JO 



— r)s 
kr 



(35) 

is essentially 



(47r)2A;2 

The reader will recognise that the argument in the exponential in equation 
the ratio of the size of the Fresnel zone to the correlation length of the fluctuations. The 
exponential takes account of the diffraction effects; its argument vanishes if the size of the 
inhomogeneities is much bigger than the Fresnel zone. The geometrical acoustics result is 
recovered by expanding the exponential in the previous equation. This is valid only if 



g^(s — r)s/kr <^ 1, 



(36) 



which is the well-known Fresnel condition. Then the integral in (|35|) is simply 



{^{) - -{2k)-'r / C(^') dr' 



(37) 

which is the geometrical acoustics result, and, unlike (|35|) , does not account for diffraction. 

Here, for simplicity, I confine myself to isotropic turbulence so that the three-dimensional 
Fourier transform becomes 



An 

C{q) = — J ((r)rsmqr dr. (38) 

This now completes the description of the propagation of waves through a homogeneous tur- 
bulent medium, the main result being equation (^). This subsection serves to illustrate the 
evaluation of equation ([ISD , which is most clearly demonstrated in a homogeneous medium. 
The Sun is, of course, inhomogeneous; hence in the following subsection I will generalise this 
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formalism to inhomogeneous media following the same principles and methods that were 
outlined above. 

2.2 Inhomogeneous medium 

For an inhomogeneous medium the correlation function C is a function of position as well 
as relative separation, i.e. ( = (^{r,TVj. Furthermore, in an inhomogeneous medium the 
amplitude is no longer simply proportional to the inverse of the distance from the source, as 
for spherical waves in a homogeneous medium. In general, one can write the Green function 
as follows: 



^(r, ri) = A{r, ri) exp[zS'(r, ri)]. 



(39) 



where the phase is given by 

S{r,ri) = J^\[r{s)]-ds . (40) 

It was shown by Munk & Zachariasen (1976), for example, that the 'normalisation' can be 
expressed as 



A(r,R) 



Air 



det 



_d d_ 

dr±idR± 



-5(r,R) 



1/2 



(41) 



where the derivatives are evaluated at r± = = 0. Now, the exponent (\E'^) reads 



(^2^ = J d^K 



A{r, R)A(R,0) 
^(r,0) 



J rf^f C(f , R) exp[i(5(r, R + f/2) + S{r, R - ?/2) 



+ S{R + f/2, 0) + S{R - f/2, 0) - 25(r, 0))]. (42) 

Expanding the exponent in powers of f , the linear terms vanish, so that one is left with the 
zero- and second-order term in f . Thus one can write 



(^?) ~ J d^R 



'A{r,R)A{R,Oy 



exp[2i(^(r,R) + 5(R,0) 



- 5(r, 0))] I d^i C(f , R) exp[l^f,r~,a,(R)], (43) 



where Cij denotes the 'phase curvature' defined as 



Cij{R) 



_d d_ 



.[5(r,R) + 5(R,0)] 



(44) 
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The argument of the exponential in (43) is summed over the indices i and j, each of which 
can take the values 1, 2 and 3. The phase curvature, Cij can be visualised as follows: Imagine 
displacing a ray which leads from to R and from R onwards to r by a small amount 5R so 
that the perturbed ray consists of a segment from to R + 6H and a segment from R + (5R 
to r. The second derivative of the phase measured at r due to a source at the origin with 
respect to perpendicular displacements at a point R is called phase curvature. Since the 
ray follows a path of stationary phase, one finds that Cij = if i 7^ j. Moreover, one can 
choose coordinates relative to the stationary ray from to r such that the only non-zero 
components of C are the second derivatives in two perpendicular directions to the stationary 
ray at any point. The phase curvature is closely related to the concept of Fresnel zones and 
will be discussed further in Seclj. 

Evaluating the integral over (i^R by the method of stationary phase, selects the unper- 
turbed path from to r as in the homogeneous case. The integral over d^r can again be 
simplified by introducing the Fourier transform of (: 



{^D ^ - r ds 
Jo 



A{r,R{s))A{Ris),Oy ' 
A(r,0) 



^-1/2^-1/2 
'-'11 '-^22 



d'q^{s)C[^^{s),R{s)]exp[tql^C-/{K{s))], (45) 

where the argument of the exponential is summed over j {j = 1,2). It is found that it is an 
excellent approximation, instead of using the normalisation of the Green function as given 
by equation (p]), to simply write A{r, R) = (47r|r — R|)^^, as in the homogeneous case. The 
deviations from homogeneity become only important in the phase. It is verified numerically 
(by using the correct Green function for a polytropic envelope, see Fan et al. 1995) that 
equation (45) can approximately be written as 

(v[/2) ^ -{An)-' J^^ ds k{s)-' J d\^{s)C[^^{s),R{s)] x exp[tql^Cr^{R{s))], 

(46) 

which is the expression that one would obtain using the amplitudes of the Green function for 
a homogeneous medium. Again equation (^6|) is much more accurate than the corresponding 
expression in the geometrical limit, which is given by 
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{^D ^- f^ds [2k{s)]-^ J df\\ C(f II, s). (47) 

Equation (45) is the most important equation of this paper. Within the assumptions made 
in its derivation, it describes the interaction of waves with the turbulent medium. My aim 
is now to evaluate (45) for acoustic waves in the Sun. In the following sections I will seek 
expressions for the various ingredients of (45), such as the scattering potential, the phase 
curvature and the correlation of the turbulence. 



3 AN INHOMOGENEOUS WAVE EQUATION 

First I derive the wave equation for solar acoustic oscillations in the presence of turbulence. 
This involves dividing the physical quantities into statistical averages and random fluctua- 
tions. The averages determine the mean properties of the envelope, whereas the fluctuations 
describe the turbulent convection and form the inhomogeneous term in the wave equation. 
For a sketch of the derivation of the linearized adiabatic wave equation in the absence of 
turbulence the reader is referred to the appendix. 

Assuming that the Mach number of the turbulent flow is small, one can linearize in 
the perturbed variables. For simplicity I neglect the advection of the waves by the turbulent 
velocity and solely consider the effect of the perturbations in the sound speed. The perturbed 
quantities are denoted by a prime. As usual, c denotes local sound speed, p matter density 
and p pressure. Thus equation (|A6|) becomes 



2^p _ Dp 



with equations ( [Al| ) and ( |A5|) unchanged. Now I repeat the steps outlined in the appendix 
to eliminate p' and p' to find an equation for the scalar x •= div^, where ^ denotes the 
displacement of the fluid. Hence, in the presence of sound speed perturbations to first order 
in c', equation (|A8| ) becomes 



= V\c\ - 90^ ■ n) - V(rx) ■ n + V\2c c'x) - V{H~^2c c'x) ■ n (49) 

Again making the substitution c^x = P~^^'^''P^ ^ind neglecting the buoyancy frequency, one 
obtains 
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+ cu^j - VV = 25c + 4(n ■ V(5c)(n ■ VV^) 

+ (2V'(<5c) - 25c cj.Vc')^, (50) 
where 5c is the fractional sound speed fluctuation c'/ c. On the Ihs I have written the familiar 
wave equation, modified by the cut-off frequency, tUc, in the absence of sound speed perturba- 
tions, and on the rhs I have written all the scattering terms proportional to the perturbation 
5c. On the rhs I can now substitute the unperturbed values for the spatial derivatives of ip by 
treating the waves as locally plane, i.e. V'^ip = —k'^ip = —'il){uj'^ — and n ■ Vip = ik^ip. 

This brings the wave equation into the form of equation (|l]) with a rhs of the form €¥{1)1^. 
The treatment of the waves as locally plane is based on a wave-like decomposition of the 
normal modes. In this approximation the modes are represented as standing waves which 
in turn are formed by mutually interfering inward and outward propagating waves. Strictly 
speaking, this approximation is only valid when the order of the mode is large and many 
wavelengths fit in any scale height of the background medium. In practice, however, it turns 
out that this is a good approximation even for moderate order. For further reference see the 
appendix of this paper and Gough (1993). 

Before proceeding to apply this formalism to the propagation of acoustic waves in the 
solar convection zone, I should recapitulate the assumptions made in this section: Plane- 
parallel geometry was assumed; as the rays that are most affected by convection are those 
confined to the surface layers of the Sun, they will feel the curvature of the Sun only as 
a small perturbation. Hence this assumption is fairly accurate. As customary, the gravi- 
tational acceleration is treated as a constant and the buoyancy frequency, A^, is ignored. 
Furthermore, perturbations in the gravitational acceleration g were neglected (" Cowling" - 
approximation). For linearization to be valid, the sound speed fluctuations and the Mach 
number of the turbulent flow are assumed to be small. Here, I only considered the scattering 
by the sound speed perturbations, i.e. by the perturbations in the refractive index. For a 
complete treatment of the interaction of waves with the turbulence one would also have 
to include those terms that arise from the advection by the turbulent velocity. Finally, for 
simplicity, I ignore the spatial derivatives of the pressure scale height H in the expression of 
the cut-off frequency and simply write it as ujc = c/2H. 
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4 FRESNEL TUBES 

In this section the phase curvature, which was introduced in Sec. 2, is calculated. At each 
point P on the wave path one can calculate how the phase of the wave due to a source at 
a fixed point Q varies when the raypath is slightly displaced from its stationary path. The 
first derivative of the phase with respect to small displacements of the ray path is zero, 
as expected by Fermat's theorem. The second derivative, also called 'phase curvature', is 



defined in equation (0). It is inversely proportional to the square of the linear extent of the 



Fresnel zone which was introduced in Sec. 2. 



Now, one can calculate the Fresnel tube for sound waves in the Sun. For simplicity one 
may assume that the solar envelope can be described by a plane-parallel polytrope. Then 
(? = CqZ is linear with depth and Cq is a constant. Indeed it is found that a polytrope of 
index /i ~ 3 yields a good fit to the observed frequencies. I should remark that this is not 
globally a good fit to the envelope, which is closer to a polytrope with fi = 3/2, but rather 
a consequence of the fact that the frequencies are dominated by the uppermost layers of the 
Sun. The ray equations for a plane-parallel polytrope can easily be solved and I briefiy quote 
the results here. If x denotes the horizontal coordinate and z depth below the surface, a ray 
that originates at the origin [x = z = 0) and propagates in the positive x-direction satisfies 



X = a 



sm-\z/aY/^ - (z/ay/^l - z/ay/^] , (51) 

where a is the depth of the lower turning point, a = u'^/c^kl, where u is the angular frequency 
of oscillation and is the horizontal wavenumber and a constant. The phase is given by 



= 2cua^/2^o ' sm-\z/ay/\ (52) 
So for the full arc between two photospheric refiections the ray has traversed a horizontal 
distance x = na and acquired a phase (p = 2tu{nxY^'^ /cq. In order to find the phase curva- 
ture at any given point along the ray, one perturbs the ray at x = Xq, while keeping the 
endpoints of the ray fixed. The second derivative of the total phase with respect to this 
displacement of the ray is the phase curvature as introduced in Sec. 2. The vertical extent of 
the first Fresnel zone is given by y/27i(d'^(f)/dz'^)~^^'^. In Fig. 3 the extent of the first Fresnel 
zone in the direction perpendicular to the ray is shown as a function of depth for a ray with 
a = O.OOli?0. The Fresnel zone increases with depth until approximately 2/3 of the depth 
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of the lower turning point from whereon it shrinks again to a smaller value at the lower 
turning point (see also Jenson et al. 1998). 

However, it should be noted that the lower turning point is a caustic of the ray where 
the asymptotic ray theory breaks down. Nevertheless, this does not considerably affect the 
calculation presented here, since the damping of the waves occurs predominantly in the up- 
permost layer of the solar envelope. In deeper regions, the turbulent Mach number is too 
small to have any effect. If one assumes that the convective cells have a typical size of a 
pressure scale height, one can note that the size of the Fresnel zone is of the same order 
of magnitude as the size of the convective cell. Therefore, diffractive effects are important 
when the acoustic waves interact with the convection. 



Since the Sun is horizontally stratified (i.e. the sound speed being a function of z only), 
the phase curvature in the direction perpendicular to the plane of the ray is the same as for 
a ray of the same length in a homogeneous medium. The phase curvature for this case was 
already calculated at the beginning of this section, and hence I can write 



2ka , . 

(53) 



z=zo s{zo)[2a - s{zo)y 
where 2a is the total length of a ray between two photospheric reflections (in a plane-parallel 
polytrope), and s{zq) is the length of the arc from to t{zo). 



5 SIMPLE MODELS OF CONVECTION 

The simplest parametrization of convection is the mixing-length theory (MLT) first devised 
by Taylor (1915, 1932). It describes turbulence by a single length scale 1^ which represents 
the dominant scale of coherent motion. The mixing length can be the characteristic height 
of a convective cell or a mean free path. One imagines either an ensemble of rising and falling 
turbulent elements that break up after having traversed a distance or a set of eddies of 
typical diameter l^. Vitense (1953) adopted a mixing length that was proportional to the 
pressure scale height, /m = aH, where H = —{dlnp/dr)~^ is the pressure scale height and 
a is the proportionality constant which has to be determined by calibrating the theory. This 
is done by evolving solar models with different values of a. The resulting one parameter 
family of mixing-length models is then calibrated against the known radius of the Sun. Thus 
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Gough & Weiss (1976) found a value of a = 1.1. 



Ignoring pressure fluctuations the convective heat flux can be written as 



F, ~ pCpwT', (54) 

where the overbar denotes horizontal average. Now consider the dynamics of the buoyant 
fluid. The vertical component of the turbulent velocity w can be estimated by equating the 
work done by the buoyancy forces as the fluid rises through a height to the kinetic energy 
gained by the fluid, i.e. 



pw^c::^\p'\gl-gpSUTyT, (55) 

where 6 = —{din p/dlnT)p. Depending on the geometry of the flow considered there may 
be another factor of order unity in equation (^5]). 

When the convection is very efficient, thermal diffusion may be ignored. The temperature 
fluctuations are then given by 



\T'\ ~ /3U (56) 

where (3 = — [dT/dz — {dT/dp)a,ddp/dz] is the superadiabatic gradient. Again I am neglecting 
factors of order unity depending on the particular theory. Combining equation (|56|) with (|55D 
gives 



w' ^ g6pfjT. (57) 

To calculate the convective heat flux one notes that the temperature fluctuation is positive 
for rising fluid and negative for falling fluid, so that all motion contributes positively to F^. 
So combining (|56D and (E^) with (|5lD one obtains 



F, ~ pCpwV = All, (58) 

where A = pCpp^/'^{gS/Ty/'^. This result was first obtained by Prandtl (1932). 

The sound speed perturbations are given by \c'\/c = \T'\/2T. \T'\/T can be calculated 
using equation (^), where the values of w and Zm are taken from a detailed model; the 
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resulting sound speed perturbations are shown in Fig. 4. 

Now one can extend the simple mixing-length approach beyond the assumption that all 
eddies are of the same size, by introducing an eddy spectrum or distribution, $(/m), that 
describes the density of eddies of size /m. In mixing-length theory $(/m) = ^(^m — ^o); where 
lo now denotes the mixing length of the standard MLT. Since there is no fundamental theory 
of turbulence which would apply to conditions inside the Sun, one has to choose a convenient 
spectrum. However, this eddy spectrum has to be normalised to yield the same convective 
heat flux as obtained by the calibrated mixing-length result. Therefore, Fc is integrated over 
the entire eddy distribution to yield a total convective heat flux of 

POO 

Fe = / Aimi^)dl^. (59) 
Jo 

For instance, one might decide to choose an eddy distribution of the form 

<I>{1^) = Bljfe-'-/'^. (60) 

Performing the integral in equation (|59|) and equating it to the calibrated MLT result, yields 
the normalisation to be i? = 1/6. 

6 RESULTS AND DISCUSSION 

Finally, I proceed to numerically evaluate the integrals in equation (45) using the phase 
curvature calculated in Sec. 4 and the turbulent spectrum given by equation (|60D . The in- 
tegrals in the correlation function e^iy {yi)V{y2)) can be simplified by expanding the terms 
that depend only on the static background state about their centre-of-mass coordinate R 
assuming that all quantities related to the static background do not vary substantially over 
a correlation length of the turbulent perturbation. 

The remaining correlation between the sound speed perturbations is then taken to be 

(5c(ri)5c(r2)) - 5c^{R)f^{r), (61) 

and similarly for the correlations involving the derivative of 5c. Here $(f) is given by equa- 
tion ( |60|) and 5c(R) is shown in Fig. 4. In general, the exponent ) is complex. The real 
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part leads to damping of the waves, and the imaginary part represents a phase shift. 



The real part of (\l'f) as a function of degree / with the lower turning point kept fixed 
is shown in Fig. 5. Since the lower turning point is kept fixed, frequency increases with 
increasing /. One can note that the real part of (^^) is negative, which means that the 
waves are damped, and the damping increases rapidly with / (or u). This does not come as a 
surprise as one would expect the scattering to become stronger with decreasing wavelength. 
One may quote, for example, scattering by small spheres (Rayleigh scattering), where the 
scattering cross section is proportional to the fourth power of k. Clearly, the scattering 
increases with the thickness of the scattering medium and thus with increasing depth of the 
lower turning point as demonstrated in Fig. 5. Shown are the results based on three different 
correlation functions for the sound speed perturbations. The solid line assumes a spectrum 
of the form given by equation (|60|). The dashed line corresponds to 

$(/^) = (2/o)-'e-'-/'«, (62) 

and the dotted line is based on 

$(/^) = /^(96/o')-ie-'-/2'«. (63) 

All three spectra are normalised according to equation (|59D . The dependence of the damping 
on the form of the convective spectrum can be relatively significant. 

Fig. 6 shows the imaginary part of (\E'^) for lower turning points of depths a = 0.01 Rq 
and a = 0.001 Rq as a function of degree. The results shown were calculated on the basis 



of the spectrum given by equation (pOl) . The phase shift is negative which implies that the 
phase of the scattered wave is advanced over the phase of the unperturbed wave. At first 
sight it may seem surprising that scattering leads to an advancement of the phase. A simple 
qualitative explanation for this result is given by Codona et al. (1985). They showed that 
a continuous random medium can cause an average advance of the arrival time of a pulse, 
and I will briefly summarise their argument in the following. 

Consider a wave propagating through a homogeneous medium from a point S to a point 
R, and assume that the random medium is concentrated in a "phase screen" at a distance s 
from the source. This screen has the effect of shifting the time of the wavefront by a small 
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random amount t{x), where x is the position on the screen and t{x) is a random distribution 
with zero mean. The travel time for a path through x is 

2 

r{x) = To + ^^-tix), (64) 

where tq is the travel time of the unperturbed wave, p = r/s{r — s), r being SR, and c the 
wave speed in the medium. By Fermat's principle the ray assumes a path such that r is 
stationary. Expanding t{x) as follows 

t{x) =to + fx + \t"x^ + (65) 

one finds that the ray crosses the screen at Xg = pet' . The travel time of the ray is then given 
by 

r(xs) ~ To + \cpt''^ - to - cpt'"^ ■ (66) 

Since t, and hence to and t', are by construction random variables with zero mean, the only 
contribution to the average travel time comes from the squared terms. The first term is 
positive corresponding to a delay, and represents the effect of geometry: The perturbed path 
is geometrically longer than the unperturbed one. The second term (called Fermat term) 
is negative and corresponds to a pulse advance. It is also twice as big as the first one re- 
sulting in an overall pulse advance. One can imagine that the rays governed by Fermat's 
principle seek out regions with a pulse advance. Thus the average travel time is given by 
r = To — The phase shifts displayed in Fig. 6 are too small (corresponding to about 

1/lOOth of a second at a frequency of 3 mHz) to be detected directly as a time delay in a 
time-distance analysis. 

Having calculated the change of amplitude of the wave between two photospheric reflec- 
tions one can convert this attenuation into a damping rate, rj, using the equation A/Aq = 
e~^'^ , where r is the travel time for a single traverse of the ray, and A/Aq the fractional 
change in amplitude per skip. Thus one finds that rj = — r^^ In A/Ao, where the amplitude 
ratio A/Aq is given by Re[exp(|(\l/f))]. Observationally this damping manifests itself as a 
line width in the measured frequency of a solar eigenmode. It is straightforward to show 
that the line width (full width at half maximum) of the acoustic eigenmodes of the Sun, F, 
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is equal to twice the damping rate rj. 

With tlie advent of long time-series of solar oscillation data as obtained, for example, 
by the ground-based GONG network and the Michelson-Doppler Interferometer on board 
the spacecraft SOHO, the line widths of high-degree modes can be measured with good 
accuracy. Recently, line widths of high-degree modes were reported by Duvall, Murawski, & 
Kosovichev (1998) and their results for the pi modes are reproduced in Fig. 7 a. In Fig. 7 
b I have plotted the fractional contribution to the line width as predicted by the method 
presented here. The results shown in Fig. 7 b are based on the spectrum given by equation 
(|6^). It is found that scattering by sound speed perturbations contributes ~ 15% to the line 
width of the high-degree modes. However, the model presented here can only be regarded 
as a toy model since the spectrum of the turbulence and the model of the surface layers 
of the Sun were chosen for convenience rather than for accuracy. But it is interesting to 
note that this model correctly predicts the variation of the line width with uo and /. The 
fractional contributions shown in Fig. 7 b are constant within the errors of the measurements. 

Another caveat that I should mention concerns the assumed structure of the turbulence. 
In Sec. 5 the convection has been described by a local mixing-length theory which presumes 
that the convection can be characterised by a locally defined spectrum. This description of 
turbulence was used in the computations in order to demonstrate the capabilities of this 
method. However, recent numerical simulations and laboratory experiments suggest that 
the local mixing-length picture is a poor description of convection in stellar envelopes (e.g. 
Nordlund & Stein 1996). Instead there are indications that the turbulent flow is extremely 
non-local, non-isotropic and hence poorly represented by a local power spectrum. Some sim- 
ulations show that heat is transported primarily through thin threads that extend all the 
way from the bottom to the top of the solar convection zone. Nevertheless, since Rytov's 
technique is an even better approximation if the perturbations are confined to thin threads, 
the method of this analysis will in principle remain valid and useful. 
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7 CONCLUSIONS 

To summarise, I have presented an analytical description of the interaction of waves with 
the turbulent medium following a similar treatment used in oceanography. The inhomoge- 
neous wave equation was solved using Rytov's technique, and several of the integrals in the 
ensuing expression could be evaluated analytically using a stationary-phase approximation; 
the remaining integrals could be solved numerically. Thus, in the framework of ray theory, 
one obtains a quantitative assessment of the scattering contribution to the line width of 
acoustic modes. 

Previous studies of the line widths of solar modes have focused on radial modes and 
generally neglected the effect of diffraction. However, diffraction becomes important for 
shallow waves for which the size of the Fresnel zone is comparable with the size of the 
convective cells. 

Scattering is not the only source of line widths. Other sources include non-adiabatic ef- 
fects associated with the radiative and convective energy transport, mechanical absorption 
of the sound waves (Gough 1980; Goldreich & Keeley 1977; Christensen-Dalsgaard, Gough, 
& Libbrecht 1989; Balmforth 1992) and partial reflection at the upper turning point (Balm- 
forth & Gough 1990). The latter showed in a simple analytical model that, for radial modes, 
partial reflection contributes significantly to the line width. 

In this paper I have isolated the effect of scattering by sound speed perturbations from 
other effects that contribute to the line width. Using a simple semi-analytical model I studied 
one of the physical processes that determine the line widths of high-degree modes. Once the 
remaining contributions to the line width are well understood chances will be excellent that 
one will be able to constrain the convective spectrum seismically. 
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APPENDIX A: LINEARIZED ADIABATIC WAVE EQUATION 

The linearized adiabatic wave equation has been derived by various authors (e.g. see Gough 
1986). In this appendix I briefly sketch its derivation, which is referred to in Sec. 3. 

Neglecting viscosity, the momentum equation for a fluid moving with velocity u can be 
written in the form 

Pj^ = - Vp + gp + ^, (Al) 
where the Lagrangian derivative D/Dt is defined as 

D d 

denotes all other forces except gravity, p is pressure and p matter density, g is the gravi- 
tational acceleration which is related to the gravitational potential (p by 

g = V0, (A3) 

where satisfies Poisson's equation 

= -AnGp, (A4) 

with G being the gravitational constant. 
The continuity equation is given by 

^ + pV-u = |^ + V-(pu) = 0. (A5) 

Moreover, I make the assumption that the changes in pressure and density are adiabatic, 
i.e. 



Dp ^ r>Dp ^ ^,Dp 

Dt p Dt Dt ^ ^ 

Here Fi denotes the first adiabatic exponent, which is defined as 

r.-f^). (A7) 



(91np^ 



s 
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and which needs to be derived from an assumed equation of state; c denotes the local sound 
speed. 



Now, equations ([A4|) , ( |A5| ) and ( |X6| ) are perturbed around the equilibrium state 

in terms of the displacement ^ of the fluid. Perturbing the equilibrium state - p = Po + p', 
P = Po + p', and = 00 + 0' ^ and by eliminating p' and p' from the momentum equation 



( |Al| ) using equations (|A4|) and (|A5|), and after some rearrangement one obtains 

= V{c\ -g^-^)- Txn (A8) 
with r = H~^(? — g, where H is the density scale height H := —{dlnp/dlnr)^^. The scalar 
X is defines as x := div^. 

Equation (|A8|) is an equation for the components of a vector. However, a vector is 
coordinate dependent so that it would be more convenient to find an equation for the scalar 
X- This is achieved by first taking the divergence of ([A8|) : 



d^X 



V^(c^X-^?n•0-n•V(^x). (A9) 



The vertical component of ^, n- ^, can be reexpressed by taking the double curl of equation 

m 

^ ■ X (V X ^) = — (n • Vx - V^n ■ |) = -gVUVx), (AlO) 



where is the horizontal Laplace operator. Solving (|A1CI|) for n ■ ^ and substituting into 
( |A9| ) yields the fourth-order equation: 

W ~ ^[^'^"'^^ " ■ ^(^"'^'^)] - N'Vlic'x) = 0, (All) 
where N is the buoyancy frequency given by 



Eliminating odd derivatives of the dependent variable by the substitution c^x = P^^^"^"^-, 
yields 

© 1999 RAS, MNRAS 000, 0-?? 



Probing Solar Convection 25 

where the quantity Uc = ^{1 — 2n- VHY^"^ is the critical cut-off frequency. For acoustical 
oscillations one can simplify equation ( |A13|) by neglecting the buoyancy frequency A^: 

(^^ + cu2^^-c2vV = 0. (A14) 

Aside from the term containing the critical cut-off frequency this is essentially a wave equa- 
tion. The cut-off frequency Uc modifies the wave equation in the sense that in the region 
where uo < the waves become evanescent. 
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Figure Al. Frcsnel zone construction. 

Figure A2. Diagram to tiie calculation of tlie Fresnel tube in tiie Sun. 

Figure A3. Extent of the Frcsnel tube (perpendicular to the ray) as a function of depth in a plane-parallel polytrope for a 
ray with a = 0.001. All lengths are measured in units of solar radii. 



Figure A4. Sound speed fluctuations in the solar envelope. 



Figure A5. Real part of —{^i) as a function of L for lower turning points of a = 0.01 Rq (a) and a = 0.001 Rq (b) assuming 
different turbulent spectra (see text). 



Figure A6. Imaginary part of as a function of L for lower turning points of a = 0.01 Rq (solid line) and a = 0.001 Rq. 



Figure A7. Panel (a) shows line widths of solar modes as measured by Duvall, Kosovichev, & Murawski (1998) using the MDI 
instrument. On panel (b) are shown the corresponding (fractional) contributions to the line width from sound speed scattering 
as predicted by my model. 



© 1999 RAS, MNRAS 000, 



p 




Figure 1: 




Figure 3: 




Figure 4: 




Figure 5: 




Figure 6: 



100 




1100 



model 

/r 



0.1 

700 



800 900 1000 1100 



Figure 7: 



